Residual Entropy of Ordinary Ice from Multicanonical Simulations 
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We introduce two simple models with nearest neighbor interactions on 3D hexagonal lattices. Each 
model allows one to calculate the residual entropy of ice I (ordinary ice) by means of multicanonical 
simulations. This gives the correction to the residual entropy derived by Linus Pauling in 1935. Our 
estimate is found to be within less than 0.1% of an analytical approximation by Nagle which is an 
improvement of Pauling's result. We pose it as a challenge to experimentalists to improve on the 
accuracy of a 1936 measurement by Giauque and Stout by about one order of magnitude, which 
would allow one to identify corrections to Pauling's value unambiguously. It is straightforward to 
transfer our methods to other crystal systems. 
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A thorough understanding of the properties of water 
has a long history and is of central importance for Ufe 
sciences. After the discovery of the hydrogen bond [1] 
it was recognized that the unusual properties of water 
and ice owe their existence to a combination of strong 
directional polar interactions and a network of specifi- 
cally arranged hydrogen bonds [2]. The Hquid phase of 
water differs from simple fluids in that there is a large 
quaUtative remnant of ice structure in the form of local 
tetrahedral ordering [3]. 

In contrast to liquid water the properties of ice are rel- 
atively well understood. Most of them have been inter- 
preted in terms of crystal structures, the forces between 
its constituent molecules, and the energy levels of the 
molecules themselves [4,5]. A two-dimensional projection 
of the hexagonal crystal structure of ordinary ice (ice I) is 
depicted in Fig. 1 (other forms of ice occur in particular 
at high pressures). Each oxygen atom is located at the 
center of a tetrahedron and straight lines (bonds) through 
the sites of the tetrahedron point towards four nearest- 
neighbor oxygen atoms. Hydrogen atoms are distributed 
according to the ice rules [2,6]: (A) There is one hy- 
drogen atom on each bond (then called hydrogen bond). 
(B) There are two hydrogen atoms near each oxygen 
atom (these three atoms constitute a water molecule). 

In our figure distances are given in units of a lat- 
tice constant a, which is chosen to be the edge length 
of the tetrahedra (this is not the conventional crystal- 
lographic definition). For each molecule shown one of 
the surface triangles of its tetrahedron is placed in the 
xy-plane. The molecules labeled by u (up) are then at 
z — l/\/24 above, and the molecules labeled by d (down) 
at z = — 1/\/24 below the xy-plane, at the centers of 
their tetrahedra. In our computer simulations informa- 
tion about the molecules will be stored in arrays of length 
N, N being the number of molecules. 

Essentially by experimental discovery, extrapolating 
low temperature calorimctric data (then available down 




FIG. 1. Lattice structure of one layer of ice 1. The up (u) 
sites are at z — l/\/24 and the down (d) sites at z = — 1/\/24. 
For each site three of its four pointers to nearest neighbor sites 
are shown. 



to about 10°K) towards zero absolute temperature, 
was found that ice has a residual entropy [7] : 



Sa = k In(VF) > 



(1) 



where W is the number of configurations for N molecules. 
Subsequently Linus Pauling [6] derived estimates oiW = 
(IFi)^ by two approximate methods, obtaining 



Pauline 



3/2 



(2) 



in each case. W — (Wi)^ is the number of Pauling 
configurations. Assuming that the H2O molecules are 
essentially intact in ice, his arguments are: 

1. A given molecule can orient itself in six ways sat- 
isfying ice rule B. Choosing the orientations of all 
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molecules at random, the chance that the adjacent 
molecules permit a given orientation is 1/4. The 
total number of configurations is thus W = (6/4)^. 

2. Ignoring condition B of the ice rules, Pauling allows 
2^^ configuration on the hydrogen bonds between 
adjacent oxygen atoms: Each hydrogen nucleus is 
given the choice of two positions, near to one of the 
two oxygen atoms. At one oxygen atom there are 
now sixteen arrangements of the four hydrogen nu- 
clei. Of those ten are ruled out by ice rule B. This 
condition for each oxygen atom permits 6/16 = 3/8 
of the configurations. Accordingly, the total num- 
ber of configurations becomes W = 2^^ (3/8)^. 

Equation (2) converts to the residual entropy 

^Pauling ^ 0.80574 .. . cal/deg/mole (3) 

where we have used R = 8.314472 (15) [J/deg/mol] for 
the gas constant as of Aug. 20, 2006 given on the NIST 

website [8] (relying on CODATA 2002 recommended val- 
ues). This in good agreement with the experimental es- 
timate 

^experimental ^ g_g2 (5) cal/dcg/molc (4) 

which was subsequently obtained by Giauque and 
Stout [9] using refined calorimetry (we give error bars 
with respect to the last digit (s) in parentheses). 

Pauling's arguments omit correlations induced by 
closed loops when one requires fulfillment of the ice rules 
for all atoms, and it was shown by Onsager and Dupuis 
[10] that Wi = 1.5 is in fact a lower bound. Onsager's 
student Nagle used a series expansion method to derive 
the estimate [11] 

lyNagle ^ 150685 (15) , (5) 

or 

^Nagie ^ 0.81480 (20) cal/deg/mole . (6) 

Here, the error bar is not statistical but reflects higher 
order corrections of the expansion, which are not entirely 
under control. The slight difference between (6) and the 
value in Nagle's paper is likely due to improvements in 
the measurements of Avogadro's number [8]. The only 
independent theoretical value appears to be one for cubic 
ice, which is obtained by numerical integration of Monte 
Carlo data [12] and in good agreement with Nagle [11]. 

Despite Nagle's high precision estimate there has ap- 
parently been almost no improvement on the accuracy of 
the experimental value (4). Some of the difficulties are 
addressed in a careful study by Haida et al. [13]. But 
their final estimate remains (4) with no reduction of the 
error bar. We noted that by treating the contributions in 
their table 3 as statistically independent quantities and 
using Gaussian error propagation (instead of adding up 



the individual error bars), the final error bar becomes 
reduced by almost a factor of two and their value would 
then read So = 0.815 (26) cal/deg/mol. Still Pauling's 
value is safely within one standard deviation. Modern 
electronic equipment should allow for a much better pre- 
cision. We think that an experimental verification of the 
difference to Pauling's estimate would be an outstanding 
confirmation of structures imposed by the ice rules. 

In this paper we provide a novel high-precision numer- 
ical estimate of 5*0 for ordinary ice. Our calculations 
are based on two simple statistical models, which re- 
flect Pauling's arguments. Each model is defined on the 
hexagonal lattice structure of Fig. 1. 

In the first model, called 6-state H2O molecule model, 
we allow for six distinct orientations of each H2O 
molecule and define its energy by 

E = -Y,h{b,slsl) . (7) 

b 

Here, the sum is over all bonds b of the lattice and {si 
and si indicate the dependence on the states of the two 
H2O molecules, which are connected by the bond) 

Hb, slsl) = ii ^ hydrogen bond, 
^ ' ''^ 1 otherwise. ^ ' 

In the second model, called 2-state H-bond model, we 

do not consider distinct orientations of the molecule, but 
allow two positions for each hydrogen nucleus on the 
bonds. The energy is defined by 

E = -J2fis^bl,bl,blbt), (9) 

s 

where the sum is over all sites (oxygen atoms) of the 
lattice. The function / is given by 

f{s,bl.bi.hl,l4)= (10) 
2 for two hydrogen nuclei close to s, 
1 for one or three hydrogen nuclei close to s, 
for zero or four hydrogen nuclei close to s. 

The groundstates of each model fulfill the ice rules. At 
P = the number of configurations is 6^ for the 6-state 
model and 2^^ for the 2-state model. Because the nor- 
malizations at /? = are known, multicanonical (MUCA) 
simulations [14] allow us in either case to estimate accu- 
rately the number of groundstate configurations [15]. Su- 
perficially both systems resemble Potts models (see, e.g., 
[16] for Potts model simulations), but their thermody- 
namic properties are entirely different. For instance, we 
do not find any sign of a disorder-order phase transition, 
which is for our purposes advantageous as the MUCA es- 
timates for the groundstate entropy become particularly 
accurate. This absence of a bulk transition does not rule 
out long-range correlations between bonds of the ground 
state configiirations, which are imposed by the conserva- 
tion of the fiow of hydrogen bonds at each molecule. In 
that sense the ground state is a critical ensemble. 
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TABLE I. Simulation data for Wi. 



N 




Tin, 


Tlz 


6-sta,te model 


2-sta,te model 


O 










• f 1 ^ eye 


* y L ' eye 




128 


4 


8 


4 
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1.51522 (49) 223 


1.51546 (15) 1096 
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12 
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1.51264(18) 503 


1.51279 (10) 1530 


0.47 


896 
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16 
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1.51075 (16) 208 


1.51092 (06) 2317 


0.32 


1600 


8 


20 


10 


1.50939 (09) 215 


1.50945 (05) 619 


0.56 


oo (fit) 


1.50741 (33) 


1.50737 (17) 


0.91 



Using periodic boundary conditions (BCs), our simu- 
lations arc based on a lattice construction set up earlier 
by one of us [17]. Following closely the method out- 
lined in chapter 3.1.1 of [16] four index pointers from 
each molecule to the array positions of its nearest neigh- 
bor molecules are constructed. The order of pointers 
one to three is indicated in Fig. 1. The fourth pointer 
is up the .2-direction for the u molecules and down the 
^-direction for the d molecules. The lattice contains 
then N = Uxnyriz molecules, where n^, Uy, and are 
the number of sites along the x, y, and z axes, respec- 
tively. The periodic BCs restrict the allowed values of Ux, 
Uy, and riz to n,j. = 1, 2, 3, . . ., Uy = 4, 8, 12, . . ., and 
Uz = 2, 4, 6, . . . . Otherwise the geometry does not close 
properly. Using the inter-site distance roo = 2.764 A 
from Ref. [5], the physical size of the box is obtained 
by putting the lattice constant a to a = 2.257 A, and 
the physical dimensions of the box are calculated to be 
Bx = 2nxa, By = {nyVi/2)a, B^ = {nz4:/VQ)a. In 
our choices of Ux, Uy, and Uz values we aimed within 
reasonable limitations at symmetrically sized boxes. 

Table I compiles our MUCA Wi estimates for the lat- 
tice sizes used. In each case a Wang-Landau recursion 
[18] was used to estimate the MUCA parameters for 
which besides a certain number of cycling events [16] a 
flatness of iJmin/^fmax > 0.5 was considered sufficient for 
stopping the recursion and starting the second part of 
the MUCA simulation with fixed weights {H[E) is the 
energy histogram, i?min is the smallest and ifmax the 
largest number of entries in the fiattened energy range). 

The statistics we used for measurements varied be- 
tween 32 X 10® sweeps for our smallest and 64 x 10^ 
sweeps for our largest lattice. Using two 2 GHz PCs 
the simulations take less than one week. The number 
of cycles, A^cyc, completed between /?i„in = and the 
groundstate are listed in the 6-state and 2-state model 
columns of table I. As each of our simulations includes 
the /3min = canonical ensemble, the (logarithmically 
coded) re- weighting procedure of chapter 5.1.5 of [16] de- 
livers estimates for Wi, which are compiled in the same 
columns. Each error bar relies on 32 jackknife bins. As 
expected, the values from both models are consistent as 
is demonstrated by Q values of Gaussian difference tests 
(see, e.g., chapter 2.1.3 of [16]) in the last column of the 
table. The 2-state H-bond model gives more accurate es- 
timates than the 6-state H2O molecule model, obviously 
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FIG. 2. Fit for Wi. 

by the reason that the cycling time, which is oc l/A^cyo is 
less for the former, because the energy range that needs 
to be covered is smaller. 

In Fig. 2 a fit for the data of the 2-state H-bond model 
to the form 



Wx{x) = Wi{Q) + axx\ x=l/N 



(11) 



is shown. The W\ = Wi (0) estimate from Fig. 2 is given 

in the the last row of the 2-state model column of table I. 
The data point for the smallest lattice is included in the 
fit, but not shown in the figure where we like to focus 
on the large region. The goodness of fit (chapter 2.8 
of [16]) is Q = 0.47 as given in the figure. Similarly the 
estimate for the 6-state H2O molecule model in the last 
row of the table is obtained with a goodness of fit Q = 
0.78. All Q values (Gaussian difference tests and fits) are 
in the range one would expect for statistically consistent 
data. The values of the fits were also consistent and 
their combined value is ^ = 0.923 (23). That we have 9 ^ 
1 reflects bond correlations in the groimdstate ensemble. 

Combining the two fit results weighted by their error 
bars leads to our final estimate 



MUCA 



= 1.50738 (16) . 



This converts into 



-.MUCA 

'0 



0.81550 (21) cal/deg/mole. 



(12) 



(13) 



for the residual entropy. 

The difference between (12) and the estimate of Na- 
gle (5) is 0.035% of the estimated Wi value (0.086% of 
Sq), which is much smaller than any foreseeable experi- 
mental error. However, within their own error bars the 
Gaussian difference test between the two estimates yields 
Q = 0.016. As the error bar in (12) covers and no sys- 
tematic errors due to finite size corrections from larger 
lattices, the small discrepancy with Nagle's result may 
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well be explained this way. In view of the large error 
bar in the experimental estimate it appears somewhat 
academic to trace the ultimate reason. 

As already (hesitatingly) pointed out by Pauling [6], 
the real entropy at zero temperature is not expected to 
agree with the residual entropy extrapolated from low 
but non-zero temperatures. In real ice one expects a 
small splitting of the energy levels of the Pauling configu- 
ration, which arc degenerate in both of our models. Once 
the thermal fluctuations become small compared with 
these energy differences, the entropy will become lower 
than the residual entropy calculated here. Such an effect 
is observed in [13] by annealing ice I at temperatures be- 
tween 85°K and 110°K. Refined models are needed to 
gain computational insights. Crossing this temperature 
range sufficiently fast allows one still to extract the resid- 
ual entropy, because the relaxation time has become so 
long that one does not have ordering of Pauling states 
during typical experimental observation times. 

It is clear that our method carries rather easily over to 
other crystal structures for which one may want to cal- 
culate residual entropies. In particular structural defects 
and impurities can be included, although one may have 
to use more realistic energy functions and lattice sizes 
put limits on low densities. Simulations very similar to 
those performed here should enable accurate estimates 
of the residual entropies for other forms of ice and var- 
ious geometrically frustrated systems [19] as well as for 
spin models in the class for which lower bounds on their 
residual entropies were derived in [20] . For more involved 
systems our approach is to design simple models, which 
share the relevant groundstate symmetries with the sys- 
tem of interest. That could, for instance, have applica- 
tions to the residual entropy of proteins by allowing for 
more realistic modeling than done in Ref. [21]. 

Finally, good modeling of water is of crucial impor- 
tance for computational progress in biophysics. Clus- 
ters of hydrogen bonds play a prominent role in water at 
room temperature. Our method allows one to calcidate 
the combinatorial factors with which small clusters 
ought to be calculated in phenomenological water models 
like those discussed in Ref. [3] . Through a better under- 
standing of hydrogen bond clusters insights derived from 
the study of ordinary ice may well be of importance for 
improving on models [22] , which have primarily been con- 
structed to reflect properties of water under room tem- 
peratures and pressures. 
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